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We  theoretically  study  the  relaxation  of  high  energy  single  particle  excitations  into  molecules  in  a 
system  of  attractive  fermions  in  an  optical  lattice,  both  in  the  superfluid  and  the  normal  phase.  In  a  system 
characterized  by  an  interaction  scale  U  and  a  tunneling  rate  t ,  we  show  that  the  relaxation  rate  scales  as 
~Ct exp[— aU2 / 1 2  In (U / 1)]  in  the  large  U/t  limit.  We  obtain  explicit  expressions  for  the  temperature  and 
density  dependent  exponent  a,  both  in  the  low  temperature  superfluid  phase  and  the  high  temperature 
phase  with  pairing  but  no  coherence  between  the  molecules.  We  find  that  the  relaxation  rate  decreases 
both  with  temperature  and  deviation  of  the  fermion  density  from  half  filling.  We  show  that  quasiparticle 
and  phase  degrees  of  freedom  are  effectively  decoupled  within  experimental  time  scales  allowing  for 
observation  of  ordered  states  even  at  high  total  energy  of  the  system. 
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Ultracold  atoms  on  optical  lattices  [1]  can  be  used  for 
simulating  strongly  interacting  quantum  many-body  sys¬ 
tems  [2-4]  with  tunable  Hamiltonian  parameters.  Although 
the  main  focus  of  cold  atom  experiments  has  been  to  obtain 
the  equilibrium  phase  diagram  of  various  models,  ultracold 
atomic  systems  also  provide  a  unique  platform  to  study  the 
intrinsic  nonequilibrium  dynamics  of  strongly  interacting 
many-body  systems.  Their  low  energy  scales  and  decou¬ 
pling  from  external  environment  lead  to  long  nonequilib¬ 
rium  time  scales  over  which  the  system  can  be  followed 
without  ultrafast  probes. 

The  issues  of  nonequilibrium  relaxation  dynamics  are 
becoming  an  important  consideration  for  the  state-of-the- 
art  cold  atom  experiments  [5-11]  as  well  as  recent  pump- 
probe  experiments  with  electron  systems  [12].  The 
tunability  of  Hamiltonian  parameters  to  access  strongly 
interacting  regimes  is  one  of  the  central  attractive  features 
of  cold  atoms.  However,  an  implied  assumption  in  con¬ 
necting  the  results  obtained  on  optical  lattices  to  the  phys¬ 
ics  of  condensed  matter  systems  is  that  the  atoms  on  the 
optical  lattice  have  achieved  thermal  equilibrium  at  low 
temperatures  after  tuning  the  parameters.  Hence,  it  is  im¬ 
portant  to  understand  the  relaxation  dynamics  and  associ¬ 
ated  equilibration  time  scales  [13-20]  of  these  systems. 

The  attractive  Hubbard  model  on  optical  lattices  is  a 
lattice  implementation  [21]  of  BCS-Bose-Einstein  con¬ 
densation  crossover  [1,22],  which  is  a  paradigm  for  under¬ 
standing  strongly  interacting  superfluids.  The  strong 
coupling  physics  is  governed  by  formation  of  tightly  bound 
molecules  which  undergo  Bose-Einstein  condensation  at 
low  temperatures.  In  this  Letter,  we  study  the  relaxation 
dynamics  of  the  attractive  Hubbard  model  in  the  strong 
coupling  limit,  where  most  fermions  are  paired  to  form 
molecules.  We  focus  on  the  decay  of  excess  unpaired 


fermions  present  in  the  system  (either  due  to  an  external 
drive  like  lattice  modulation  or  due  to  sweeping  of  the 
Hamiltonian  parameters)  to  form  molecules.  For  these  high 
energy  excitations,  energy  conservation  requirements  lead 
to  a  very  slow  decay  rate  that  scales  superexponentially 
with  the  ratio  of  the  interaction  strength  to  the  bandwidth 
of  the  system.  Thus,  at  strong  coupling,  the  fermionic 
quasiparticles  and  the  motion  of  the  bosonic  molecules 
effectively  decouple  in  the  nonequilibrium  dynamics  of 
the  system.  This  leads  to  the  possibility  of  observing 
ordered  states  in  these  systems  even  at  a  high  total  energy 
of  the  system.  Using  a  particle-hole  transform  to  map  this 
problem  to  that  of  spin  mediated  decay  of  double  occu¬ 
pancies  in  the  repulsive  Hubbard  model,  we  compute  the 
decay  rate  both  in  the  low  temperature  superfluid  phase  and 
in  the  high  temperature  paired  phase  for  arbitrary  filling 
fractions  in  the  lattice.  We  find  that  the  decay  rate  de¬ 
creases  both  with  temperature  and  with  the  deviation  of 
the  fermion  density  from  half  filling  on  either  side.  We 
discuss  the  implications  of  these  results  for  maintaining 
adiabaticity  during  a  sweep  of  Hamiltonian  parameters. 

We  consider  the  one  band  attractive  Hubbard  model  for 
fermions  on  a  3D  cubic  optical  lattice 

H  =  ~  U^nnna,  (1) 

(ij)  1 

where  t  is  the  tunneling  matrix  and  U  is  the  local  attraction 
between  the  fermions.  For  large  U/t,  the  fermions  are 
paired  to  form  molecules  with  a  large  binding  energy 
—  [/  which  undergo  Bose  condensation  at  a  temperature 
=  4 t2/U,  controlled  by  the  kinetic  energy  scale  of  the 
molecules.  We  consider  two  temperature  regimes:  (a) 
the  low  temperature  (T  «  J)  superfluid  phase,  where 
the  molecules  are  Bose  condensed,  and  (b)  the  high 
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temperature  (T  ~  t  »  J)  phase  where  the  molecules  do 
not  have  phase  coherence,  and  calculate  the  decay  rate  of 
unpaired  fermions  in  both  phases. 

The  formation  of  a  molecule  from  two  unpaired  fermi¬ 
ons  is  forbidden  unless  the  binding  energy  of  the  molecule 
( —  C/)  is  carried  off  by  other  excitations  in  the  system. 
There  are  two  different  modes  of  excitations  where  this 
excess  energy  can  be  dumped:  (a)  kinetic  energy  of  other 
unpaired  fermions  (with  scale  ~t)  and  (b)  kinetic  energy  of 
the  molecules  (with  scale  —/).  In  this  Letter,  we  assume 
(b)  is  dominant,  which  limits  the  density  of  unpaired 
fermions  pQX  <  J/t  ~  t/U.  To  absorb  the  binding  energy 
ft  ~  U/J  ~  U2/t2  molecular  excitations  have  to  be  cre¬ 
ated.  The  corresponding  matrix  element  is 


where,  to  leading  order,  we  have  used  n\  ~  nn  for  large  n 
and  nJ  =  U  to  write  the  final  form.  Here  a  is  a  temperature 
and  density  dependent  constant  which  depends  on  the 
details  of  the  process.  Within  Fermi’s  golden  rule,  the 
decay  rate  T  ~  M 2  and  thus  decreases  superexponentially 
with  U/t. 

To  obtain  a  physical  picture  of  the  decay  process,  it  is 
instructive  to  use  a  particle-hole  transformation  [23]  which 
maps  the  attractive  (negative  U)  Hubbard  model  to  a 
repulsive  (positive  U)  Hubbard  model.  The  attractive 
model  with  zero  magnetization  at  any  density  is  equivalent 
to  the  repulsive  model  at  half  filling  (one  particle  per  site) 
with  a  finite  magnetization  proportional  to  the  deviation  of 
the  fermion  density  in  the  attractive  model  from  half  fill¬ 
ing,  i.e.,  m  =  (1/2)(1  —  p),  where  p  is  the  fermion  density 
in  the  attractive  Hubbard  model.  Under  this  transforma¬ 
tion,  molecule  formation  is  mapped  to  the  formation  of  a 
Mott  insulator,  and  the  unpaired  fermions  are  equivalent  to 
the  high  energy  double  occupancy  (doublon)  hole  excita¬ 
tion,  with  the  binding  energy  of  the  molecules  playing  the 
role  of  the  Mott  gap.  At  half  filling,  the  low  energy  physics 
of  the  repulsive  Hubbard  model  reduces  to  an  antiferro¬ 
magnetic  Heisenberg  model  which  exhibits  a  canted  anti¬ 
ferromagnetic  (CAFM)  order  in  its  ground  state  in  the 
presence  of  finite  magnetization.  This  spin  ordering  is 
equivalent  to  the  emergence  of  superfluidity  in  the  attrac¬ 
tive  model  with  the  kinetic  energy  of  the  molecules  playing 
the  role  of  spin  wave  fluctuations.  The  equivalence  of  the 
relevant  quantities  under  the  mapping  is  shown  in  Table.  I. 
Here,  we  will  use  the  language  of  the  spin  model  to  look  at 
quantitative  estimates  of  the  decay  time  scales. 

Decay  in  the  superfluid  phase . — The  superfluid  phase  of 
the  attractive  fermions  is  represented  by  the  CAFM  phase 
for  the  spins  in  the  repulsive  model.  As  a  hole  (doublon) 
hops  in  the  background  of  a  Mott  insulator  with  CAFM 
ordering  [shown  in  Fig.  1(a)]  from  a  site  i  to  a  site  j ,  it 
disrupts  the  spin  texture  and  creates  purely  ferromagnetic 
bonds  between  nearest  neighbors,  each  of  which  gains  an 


TABLE  I.  Equivalence  of  different  quantities  under  the  map¬ 
ping  between  the  attractive  and  the  repulsive  Hubbard  model. 


Attractive  model 

Repulsive  model 

Unpaired  fermions 

Doublon  hole  pairs 

Binding  energy 

Mott  gap 

Deviation  from  half  filling 

Magnetization 

Superfluid  order 

Canted  antiferromagnetic  order 

Kinetic  energy  of  molecules 

Superexchange  energy 

energy  of  (Jx/ 2),  where  x  =  1  —  4m2  is  proportional  to  the 
antiferromagnetic  component  of  the  spin  order.  Hopping  of 
the  hole  along  a  path  creates  ferromagnetic  bonds  in  the 
directions  transverse  to  this  path,  thus  creating  a  domain 
wall  in  the  system,  as  shown  in  Figs.  1(b)  and  1(c).  In  a 
cubic  lattice  each  hop  creates  z  —  2  =  4  broken  bonds.  In 
order  to  absorb  the  energy  [/,  the  hole  needs  to  traverse  a 
path  of  length  n  =  2 U /[(z  ~  2)Jx\.  Within  Fermi’s  golden 
rule,  the  decay  rate  is 

r  =  2Trpexv(n) \Mfi{n) \2,  (2) 

where  pQX  is  the  density  of  holes  (or  equivalently  of  un¬ 
paired  fermion  excitations),  and  Mj-fn)  is  the  matrix  ele¬ 
ment  connecting  the  initial  state  with  a  doublon  and  a  hole 
to  the  final  state  with  a  domain  wall  of  length  n 

Here,  v(n)  is  proportional  to  the  number  of  self-avoiding 
paths  of  length  n  connecting  the  doublon  and  the  hole 

/2ft 

d3r—  -  2)/xV»,  r)G(r),  (4) 

where  S{n,  r )  is  the  number  of  self-avoiding  paths  of  length 
ft  connecting  two  points  at  a  distance  r,  and  G(r)  is 
the  dimensionless  doublon-hole-pair  correlation  function. 
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FIG.  1  (color  online).  Hopping  of  a  hole  in  a  canted  antiferro¬ 
magnetic  background,  (a)  A  single  hop  moves  back  one 
spin,  creating  broken  bonds  shown  by  dashed  lines,  (b), 
(c)  Configurations  before  and  after  multiple  hops  of  a  hole. 
The  solid  line  denotes  the  trajectory  of  the  hole,  while  the  dashed 
lines  in  (c)  shows  the  broken  bonds. 
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We  assume  that  the  doublons  and  holes  are  uncorrelated, 
i.e.,  G(r)  =  1.  This  assumption  is  valid  in  the  limit  of  low 
density  of  doublons,  precisely  the  limit  we  are  interested 
in.  Then,  v(n)  =  (z-2)jxS(n)>  where  S(n ),  the  total  number 
of  self-avoiding  paths  of  size  n.  S(n )  ~  gnnk  for  large  n , 
and  the  constants  g  =  4.68  [24]  and  k  =  1/6  [25]  for  a 
cubic  lattice  have  previously  been  computed  in  the  context 
of  polymer  physics.  Using  these,  we  obtain  the  relaxation 
rate  of  doublons  in  the  repulsive  Hubbard  model  or  equiv¬ 
alently  the  relaxation  rate  of  the  unpaired  fermions  in  the 
attractive  Hubbard  model, 

Note  that,  for  the  attractive  case,  v  =  2p  —  p2  is  a  measure 
of  the  filling  factor  which  vanishes  both  at  p  =  0  (empty 
band)  and  p  =  2  (completely  filled  band)  and  attains  its 
maximum  value  at  half  filling  (p  =  1).  The  low  tempera¬ 
ture  decay  rate,  plotted  for  different  densities  in  Fig.  2(a), 
shows  the  expected  superexponential  scaling  with  U/t.  As 
we  move  away  from  half  filling,  the  energy  lost  in  a  single 
hop  decreases,  and  longer  domain  walls  are  required  to 
absorb  the  excess  energy,  leading  to  a  slower  decay  rate. 

Decay  in  the  high  temperature  normal  state. — We  now 
consider  the  decay  of  the  single  particle  excitations  in  the 
high  temperature  phase  (T  ~  t  »  /),  where  the  fermions 
are  still  paired  into  molecules  but  there  is  no  superfluidity. 


FIG.  2  (color  online).  The  decay  rate  of  unpaired  fermions 
(a)  in  the  superfluid  phase  and  (b)  in  the  high  temperature  normal 
state  as  a  function  of  U/t  for  different  fermion  densities:  (i)  half 
filling  (solid  black  line),  (ii)  pf  =  0.7  (red  dashed  line),  and 
(iii)  pf  =  0.5  (blue  dotted  line). 


In  terms  of  the  spin  model,  this  regime  corresponds  to  a 
completely  spin  disordered  phase  with  no  spatial  or  dy¬ 
namic  correlations.  The  single  site  Hilbert  space  can  be 
occupied  by  an  \  spin  with  probability  1/2  +  m  or  by  a 
down  spin  with  probability  1/2  —  m  (we  are  working  at 
fixed  magnetization). 

As  in  the  low  temperature  phase,  the  motion  of  holes 
(doublons)  pushes  the  spins  along  the  trajectory  by  one 
site.  However,  the  energy  lost  in  a  given  hop  now  depends 
on  the  configuration  of  the  neighboring  spins  along  the 
path.  Since  there  is  no  spatial  correlation  between  the 
spins,  the  energy  lost  in  each  hop  can  be  treated  as  an 
independent  random  variable  which  takes  the  values 
Jr/2  (r  =  — 4,  —  3, . . .,  4),  with  the  probability  P(r)  = 
(< a/b)r/ 2  Xto4  C,-4CI+ra'64_',  where  nCk  are  the  binomial 
coefficients,  a  =  (1/2  —  m)2,  and  ft  =  (1/2  +  m)2.  The 
mean  energy  lost  in  each  hop  is  0,  while  the  variance  of 
the  distribution  is  given  by  J2x/ 2.  By  central  limit  theo¬ 
rem,  the  total  energy  lost  in  l  steps  is  a  Gaussian  random 
variable  with  zero  mean  and  a  variance  lJ2x/ 2, 

P(E,  l )  =  J—  exp[-£2/ (J2lx)].  (5) 

\J7TXU 

Since  the  doublon  needs  to  lose  an  energy  U  to  decay,  one 
must  now  average  over  decay  processes  from  paths  of 
length  l  >  n  with  the  probability  distribution  P(U,  /). 
The  square  of  the  matrix  element  for  a  process  involving 
l  hops  is  given  by  M2~t2(t2/e2^(t2/ety---(t2/e2)~ 
t2(2t2/lJ2x)1,  where  et  is  the  energy  lost  after  i  steps  and, 
to  leading  order,  we  have  replaced  e2  by  its  average  value 
iJ2x/ 2.  To  see  the  scaling  in  t/U,  we  note  that  the 
Gaussian  probability  distribution  P(U,  l )  has  a  width 
l  ~  2U2/(xJ2),  and  hence  the  square  of  the  matrix 
element  scales  as  t2(t/U)21.  Then,  summing  over  all  paths 
with  l>  n, 

/oo  /  t  \2l 

dlgllkP(U,  /)(— j  .  (6) 

The  high  temperature  decay  rate,  plotted  as  a  function  of 
U/t  in  Fig.  2(b),  is  orders  of  magnitude  smaller  than  the 
low  temperature  decay  rate.  To  understand  this,  note  that 
the  motion  of  the  high  energy  pair  can  both  excite  and 
deexcite  low  energy  modes  in  the  background  system.  At 
zero  temperature,  these  low  energy  modes  are  unoccupied 
and  the  motion  of  the  high  energy  pair  then  leads  to 
excitation  of  these  modes.  As  temperature  increases,  the 
occupation  probability  of  the  low  energy  modes  increases 
and  the  motion  of  the  high  energy  pair  randomly  leads  to 
excitation  and  deexcitation  of  these  modes.  Thus  the  high 
energy  pair  loses  its  energy  more  efficiently  at  lower 
temperatures  and  decays  faster. 

Adiabaticity  and  sweep  rates. — In  cold  atom  experi¬ 
ments  the  strongly  interacting  regime  of  model  systems 
is  accessed  by  tuning  the  Hamiltonian  parameters  at  a  finite 
rate,  which  is  limited  by  the  lifetime  of  the  atoms  in  the 
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trap.  The  tuning  process  should  be  adiabatic  to  remain  in 
the  interesting  low  temperature  regime  for  the  system. 
Since  the  microscopic  relaxation  processes  determine  the 
time  scale  for  equilibration,  the  relaxation  time  scales, 
along  with  experimental  sweep  rates,  would  determine 
the  limits  of  adiabaticity  in  these  experiments.  We  now 
make  these  ideas  more  precise  by  looking  at  the  constraints 
due  to  the  slow  decay  of  unpaired  fermions. 

In  the  large  U/t  limit  of  the  attractive  Hubbard  model, 
the  unpaired  fermion  density  in  equilibrium  pex  ~ 
exp  {—U/T),  where  T  is  the  temperature  of  the  system. 
We  assume  an  adiabatic  sweep  of  U/t  at  a  constant  rate 
y  =  (U/t)  and  try  to  assess  the  limits  where  adiabaticity 
fails.  At  low  temperatures,  the  entropy  mainly  comes  from 
the  kinetic  motion  of  the  molecules;  so  for  a  constant 
entropy  process,  T/J  =  X/ 4  or  U/T  =  U2/(\t2)  along 
the  sweep,  where  A  is  a  constant.  Now,  adiabaticity  will 
be  maintained  in  the  regime  where 

Pex  =  ~(U/T)pex  =  Pex  ^  -r(t//f)pex-  (7) 

A  t 

As  the  microscopic  rate  T  goes  down  superexponentially 
with  U/t,  this  criterion  would  set  an  upper  limit  of 
(U  /  t)max(y),  which  is  the  maximum  U/t  up  to  which  the 
system  remains  adiabatic  with  a  sweep  rate  y.  Our  analysis 
shows  that  it  is  extremely  difficult  to  keep  the  system  fully 
adiabatic  in  the  strong  coupling  limit  when  either  the 
tunneling  or  the  interaction  is  being  changed,  as  the  re¬ 
laxation  time  scale  of  unpaired  atoms  (or  doublons  for  the 
repulsive  case)  can  be  anomalously  long.  Experimentally 
this  long  time  scale  should  manifest  itself  as  a  saturation  in 
the  molecular  fraction  with  the  saturation  occurring  at 
smaller  values  of  U/t  for  faster  sweep  rates. 

At  the  same  time,  if  relaxation  of  unpaired  fermions  is 
very  slow  (longer  then  the  time  scale  of  experimental 
measurements),  then  they  can  be  considered  as  infinitely 
long  lived  and  completely  decoupled  from  other  degrees  of 
freedom  in  the  system  like  the  phase  fluctuations  of  the 
superfluid  order  parameter.  Similarly,  in  the  repulsive 
Hubbard  model,  if  the  goal  is  to  observe  antiferromagne¬ 
tism,  one  may  worry  that  a  small  number  of  doublons  can 
release  enough  energy  to  destroy  magnetic  order.  If  the 
doublons  are  very  long  lived,  there  will  be  a  long  time  scale 
over  which  one  can  neglect  relaxation  of  doublons  and 
analyze  the  quasiequilibrium  with  “unbreakable”  dou¬ 
blons.  Thus,  within  experimental  time  scales,  there  is  an 
effective  spin-charge  decoupling  which  makes  it  easier  to 
observe  spin  ordering  even  in  the  presence  of  high  energy 
charge  excitations.  The  idea  of  realizing  metastable  states 
with  long-lived  doublons  has  also  been  discussed  in  the 
context  of  the  r\ -paired  state  in  the  repulsive  Hubbard 
model  [26]. 

We  have  studied  the  decay  of  unpaired  fermions  in  an 
attractive  Hubbard  model.  We  have  shown  that  the  decay 
rate  scales  as  ~C7exp[—  aU2 / 12  \n(U / 1)\  for  large  U/t 
and  computed  the  exponent  a  both  at  low  temperatures 


(superfluid  phase)  and  high  temperatures  (normal  state  of 
molecules).  We  find  the  decay  rate  decreases  with  increase 
in  both  temperature  and  the  deviation  of  the  fermion  den¬ 
sity  from  half  filling.  We  also  discussed  implications  of  our 
analysis  for  realizing  many-body  states  in  optical  lattices. 
The  downside  of  the  long  relaxation  times  is  that  it  is 
difficult  to  change  parameters  of  the  system  fully  adiabati- 
cally.  The  upside  of  slow  relaxation  is  that  there  is  effective 
decoupling  of  different  degrees  of  freedom.  So,  for  ex¬ 
ample,  one  may  be  able  to  achieve  equilibration  of  phase 
(magnetic)  degrees  of  freedom,  even  when  there  is  a  finite 
density  of  unrelaxed  single  fermions  (doublons). 
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